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ABSTRACT 

A series of two dimensional numerical simulations of explosive nuclear burning is 
presented for white dwarfs near the Chandraskhar mass. We assume that the burning 
begins as a slow deflagration front at or near the center of the star, and continues until 
the density in the burning regions has declined to about 10 7 g cm -3 , where the flame is 
essentially extinguished. We employ a novel numerical representation of the turbulent 
flame brush based upon ideas previously developed for modelling laboratory combustion 
and explore in some detail the sensitivity of the outcome to the manner in which burning 
is initiated. In particular, we simulate 1) a centrally ignited deflagration, 2) off-center 
ignition at a single "point" , and 3) simultaneous off-center ignition at five "points" . We 
find that the amount of 56 Ni that is produced and other observable properties depend 
sensitively upon how the fuel is ignited. Because of the immediate onset of buoyant 
acceleration, the burning region in models ignited off center rises toward the surface more 
quickly than in the (commonly assumed) case of central ignition. With the exception 
of the model that ignited off-center at a single point, all models are unbound at the 
end of the computations and between 0.59 M (central ignition) and 0.65 M (ignition 
at multiple "points") of matter are processed into nuclear burning products. These 
results would guarantee an observable, though weak, Type la supernova. Our results 
are expected to change for simulations in three dimensions, especially for the off-center 
ignitions discussed in this paper, and late detonations driven by pulsations are not 
unambigously excluded. We can, however, state that the chances for a direct transition 
to a detonation appear small because, in all our models, the turbulent velocity of the 
burning front remains very subsonic. 
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1. Introduction 

The mechanism whereby an accreting carbon-oxygen 
white dwarf explodes as a Type la supernova con- 
tinues to be uncertain. While considerable recent 
attention has been given to to the so called l: sub- 



nited flame fronts with perturbed surfaces that first 
go through the linear Raleigh- Taylor (RT) instability 
phase, and later develop rising bubbles of hot, burned 
material and turbulence ( Muller fc Arnett 1986), Livncj 
1993| , [Khokhlov 1995| , |Nicmeyer fc Hillebrandt 1995| ). 



Chandrasckhar mass models" (e.g., [ben fc Tutukov 



1991. 


Limongi & Tornambc 1991 


, Kenyon et al. 1992 , 


Munari & Rcnzini 1992 




Woosley & Weaver 1994, 



Livnc & Arnett 1996), there remain compelling rea- 
sons to think that a portion, perhaps most Type la 
supernovae are the explosion of dwarfs that have ap- 
prcjached the Chandrasckhar mass, M c h an ~ Mq. 
These reasons include: 1) the recent identification 



This assumption of central ignition simplifies the cal- 
culation and seemed reasonable based upon existing 
one dimensional models of the pre-explosive evolu- 
tion. Unfortunately, these models all share a com- 
mon fate - an inadequate amount of burning occurs 
to produce a credible explosion. This has led some to 
speculate (|Khokhlov 1995|,|Arnett fc Livne 19944 [Ar ~ 



of supersoft x-ray sources (Rappaport, DiStefano, & 
smith 1994) as the possible progenitors; 2) the re- 
quirement of conditions that exist only in this class 
of model for the nucleosynthesis of the isotopes 48 Ca, 
54 Cr, and 50 Ti flWoosley fc Weaver 1986], |Woosley fc 



Eastman 1995); 3) the good agreement of the pro- 



totypical M c han model, W7 (Nomoto, Thielemann 
k yokoi 19"f|) wiLh Ih 



bseived speclia of typi- 
cal ' Type la supernovae, 4) Ihe good agieemeiil wilh 

light curves from this class of model with a variety of 



Type la supernovae (Hoflich & Khokhlov 1996); and 



5) the difficulty achieving a robu st explosion in three 
dimensional sub-M cna n models (Benz 1996, Garcia 



5enz, Bravo, fc Woosley 1996 ). 

Despite this evidence, there remains great uncer- 
tainty both as to whether growth of the CO dwarf 
to Mchan is frequently realized and the nature of the 
nuclear explosion once it has begun (Niemeyer & 
Woosley 1996] ) . The similarity of the progenitor star 
mass, and therefore of the initial density structure, 
suggests a common outcome from this kind of event. 
However, there are properties of the progenitor star - 
the accretion rate, the original mass of the CO white 
dwarf before accretion began, its carbon to oxygen 
ratio, the metallicity, and perhaps its rotation rate - 
that might lead to diversity. Garcia-Senz & Woosley 



199, 



few); 



init ial differ* 



111 pailiculai, have suggested llial small 
elites iniglil be amplified because llie eun- 
vedlive imiaway thai pieeeeds Hie explosion may lead 

to ignition at multiple and unpredictable locations. If 
the global characteristics of the event, i.e. the ejecta 
composition and energy, depended strongly on the ig- 
nition conditions some spread of outcomes would fol- 
low as a natural consequence. 

Previous multi-dimensional simulations of M c han ex- 
plosions have studied the dynamics of centrally ig- 



nctt fc Livnc 1994b| , |Niemeyer fc Woosley 1996| ) that 



the explosion must follow one or more pulsations of 
the white dwarf in which energy is shared between 
ashes and fuel and a detonation occurs. 

Here, we take a different approach, returning to the 
simple (non-pulsing) deflagration models and consid- 
ering the effects of off center and multiple point ig- 
nition. While our immediate motivation is the work 



of GW, there are other lines of reasoning (e.g., Iben 



1982| ) that also suggest off-center ignition. In the GW 



model ignition points are set by conditions inside the 
star when its central temperature reaches 6 x 10 s K. 
At this point the time for nuclear burning and the 
convective cycle of a burning bubble become compa- 
rable. For certain circumstances, which the star may 
naturally realize as it continuously increases its burn- 
ing rate, rising bubbles of carbon continue to burn 
and become more buoyant as they float in an evo- 
lution that is distinctly non-adiabatic. GW's model 
suggests that the flame is born on the surfaces of ris- 
ing buoyant bubbles with a typical diameter of about 
10 km, at a radial distance between 100 and 200 km 
off-center. 

To follow these rising burning bubbles, we employ 
a novel prescription for flame tracking. Usually the 
macroscopic burning front is numerically represented 



speed (Livne 1993 


. Khokhlov 


Khokhlov 1995, Niemeyer &j 



Hillebrandt 1995 ). Here we use a "turbulent flame 
brush" whose properties are discussed in section (|J). 
Because the prescription is new we calculated both 
centrally and non-centrally ignitied flames for com- 
parison. 

We find that the amount of 56 Ni that is produced, 
the overall energy released, and especially the ra- 
dial distribution of the burning products are all quite 
sensitive to how ignition is simulated. However, we 
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are unable to conclude unambiguously from the two- 
dimensional studies whether the star explodes on the 
first try or pulses. Continued studies in three dime- 
sions that calculate the evolution for a longer time 
as well as a better understanding of turbulent flame 
physics will be necessary to resolve this point. 

2. Numerical representation of the turbulent 
flame brush 

Owing to the enormous difference between the 
grid resolution of multidimensional SN la simula- 
tions, A w 10 km, and the thickness of laminar ther- 
monuclear flames, /th ~ 10~ 4 cm, it is impossible 
to use the nuclear rate equations, evaluated at the 
grid averaged temperature Ta and density pa, to di- 
rectly compute the energy generation rate S of a zone 
that contains fragments of the turbulent flame brush 



("flamelets") 


(Livnc 1993, b 


Chokhlov 1993, 


Arnett & 


Livne 1994b, 


Khokhlov 1995 


, Niemeyer & Hillebrandt 



1995). In other words, one has to suitably transform 
the strongly peaked, unresolved reaction structure of 
flamelets immersed in unburned fuel into an average 
value for S of the whole zone. Furthermore, the prop- 
agation speed and thickness of the burning front have 
to be modeled in a realistic fashion. Previously, all 
this has been done by decoupling the numerical flame 
front from the grid averaged thermodynamic quanti- 
ties and using an artificial scheme to propagate a thin 
interface separating fuel and ashes with a prescribed 
"turbulent" flame speed u(A). The latter can be ob- 
tained by invoking isotropic turbulence, in which case 
a thick "turbulent flame brush" develops that prop- 
agates rough ly at the speed of the largest unresolved 
eddies v(A) ( Niemeyer fc Hillebrandt 1995 , Niemeyer 



k Woosley 19961 ). 

However, total decoupling from the thermody- 
namic state of the numerical cell prevents fluid expan- 
sion or compression from directly affecting the flame. 
Only by introducing further artificial rules can local 
ignition or quenching of the flame by temperature 
changes be included in these schemes. An alterna- 
tive way to treat the prohlcm is to conple the grid 



tcrrlpcraturc T A to the energy generation rate via, the 



burning time Tburn = tnuc/S. Obviously, if the physi- 
cal burning times were known in each zone, the code 
would be able to decide whether to assume the pres- 
ence of highly localized flamelets in the zone, or else 
to burn slowly and homogeneously as it does prior to 
flame ignition. We therefore state the main assump- 



tion of our flame model: 

^burn^ 



7b urn 

(pA,f A ,X iA ) . (1) 



This assumption is reasonably well satisfied. Immedi- 
ately after the flame front enters a grid zone, nuclear 
energy is deposited into the zone in the form of inter- 
nal energy €ia- By virtue of the small heat capacity 
c p of degenerate matter, reflected by the numerical 
equation of state, and the relation dT ~ c~ 1 dej, small 
growth of £iA already leads to a strong rise of the 
zone temperature Ta- Consequently, Ta quickly ap- 
proaches the final temperature Tb of the burned mix- 
ture, even if only a small part of the zone is burned. 
The numerical burning time as computed from the 
grid averaged temperature (the rhs of equation |l|) is 
thus mostly evaluated at Ta ~ Tb. This is consistent 
with the fact that the physical, grid averaged burning 
time (the Ihs of equation [l]) is also dominated by burn- 
ing at Tb due to the strong temperature dependence 
of S. 

In order to distinguish between low temperature, 
volume filling burning that is fully resolved by the 
grid and the high temperature regime that gives rise 
to flame formation, Tburn must be related to the char- 
acteristic length scale of burning regions. To this end, 
we invoke the stationarity assumption Tburn = r trans, 
where T tra ns stands for th e heat transport time scale 
( Landau fc Lifshitz 1991 ). Two different transport 
mechanisms compete on all scales: thermal conduc- 
tion and turbulent transport. Both can be described 
by transport coefficients that connect T tra ns (and, 
hence, Tburn) to the burning length /bum via 



P°P ^burn 



/r 



Ttu 



(2) 



The turbulent transport coefficient is the so-called 
"eddy viscosity" utm oc v(A) A, where v(A) again de- 
notes the mean magnitude of turbulent velocity fluc- 
tuations on the grid scale. One possible technique 
to compute v(A) and z^tur with the help of a subgrid 
(SG) model for the unresolved tur bulent kinetic en- 
ergy q g v(A ) 2 /2 was employed by Niemeyer fc Hille- 



brandt 1995 for SN la calculations; the same method 
was applied in this work. One must be cautious, how- 
ever, to evaluate the eddy viscosity at the burning 
length scale, that is 



^tur bin 
Ktur(A) 



burn 



4/3 



(3) 
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where we have used the Kolmogorov scaling law 
v(l) oc I 1 / 3 . Turbulent transport is realized on larger 
scales while at the same time, thermal conduction 
dominates small scale diffusion owing to the differ- 
ent scaling behaviour of both transport coefficients. 
Whereas vtm depends on turbulent velocity fluctua- 
tions that decay on small scales, thermal conduction 
by electrons takes over on microscopic scales. The 
position of this transition is defined by the burning 
time (and, therefore, by temperature). Equations (^) 
and (g) yield the burning lengths at given Tbum cor- 
responding to both transport mechanisms: 



n3/2 A _ a 



, . .1 ,, 1 — , Ztur — (Tburn' y tur(A))' :>/ A 

Cpp 



(4) 

Fig. demonstrates the behaviour of both quan- 
tities at a typical turbulent intensity and density for 
thermonuclear flames in white dwarfs. Since the mi- 
croscopic burning time at Tb « 8 x 10 9 K in the 
flamelet regime is T burn ~ 10~ 12 s, the graph shows 
that microscopic flame dynamics are governed by 
thermal conduction. 

Phcnomcnologically, the burning regime that is re- 
alized at a given temperature is the one with the max- 
imum burning length: 



max(Z cond ,Z t ur) 



(5) 



If ^burn ^> A, the burning region is fully resolved and 
the energy generation rate is given directly by the nu- 
clear reaction rate. This only occurs prior to ignition 
of the zone by the spreading flame brush. In case 
of I 

bum < A, the energy generation rate can only 
be computed by means of the effective flame speed 
it(A), for reasons explained above. Two further sub- 
cases are possible: if Ibaxn — hm, the zone is not yet 
hot enough to burn in the flamelet regime, as tur- 
bulent transport still dominates at the burning scale. 
We can approximate the effective "flame speed" in 
this transition regime (corresponding to the preheat- 
ing region of laminar flames) by u(A) « Zburn/ T bum- 
As soon as the flame brush fully enters the zone, the 
temperature rises to the point where flame transport 
is governed by thermal conduction on the scales de- 
fined by the burning time, i.e., ibum = ^cond "C A. 
The effective macroscopic burning speed then follows 
from the maximum of all competing transport veloc- 
ities at the grid scale: 



where U\ am represents the laminar flame speed, adopted 
from Timmcs & Woosley 1992. Here, u r t(A) is the 
asymptotic velocity of rising buoyant bubbles with 
the radius A: 



«rt(A) =By/g cS A 



(7) 



( Layzer 1955 ), where B = 0.511 is a constant, g e g = 
At g is the effective gravitational acceleration, and the 
Atwood number At — (p u — pb)/(pn + pb) is a mea- 
sure for the density contrast of burned (b) and un- 
burned (u) material. Khokhlov 1995| confirmed that 
RT-unstable flames in tubes with the radius A propa- 
gate with a speed given by equation (0) with the help 
of three-dimensional simulations. In our calculations, 
however, we used a smaller coefficient (B = 0.25) than 
the single-bubble value derived by Layzer in order to 
account for the presence of many bubbles with smaller 
radii in one cell. The turbulent fluctuation velocity 
u(A) is obtained directly from the kinetic SG-energy: 



«(A) = ^JYq 



(8) 



The effective energy generation rate is then calculated 
from 



1 



pS )a = as; d xpS — a — ■ (9) 

Assumption (|l|) enables us to determine the pres- 
ence of flame segments in a given zone. Therefore, no 
artificial prescription for the thickness of the flame 
brush is required, since its structure is now fully de- 
termined by the temperature and fuel mass fraction 
of the fluid. Moreover, the flame is now allowed 
to spontaneously ignite if the temperature increases 
sufficiently (e.g., by compression) to yield a burning 
time that indicates the formation of localized energy 
sources. 

In addition to the energy generation rate the model 
has to provide proper propagation of the flame sur- 
face. A possible method follows in a very natural way 
from the fact that temperature is used as an indica- 
tor of flame presence. The diffusion of internal energy 
by unresolved turbulent motions is equivalent to ther- 
mal conduction with a modified transport coefficient 
which is of the same order of magnitude as the eddy 
viscosity. Since the eddy viscosity is known from the 
SG-model, turbulent energy diffusion can be param- 
eterized by 



u(A) = max(M lam ,t;(A), v rt (A)) , 



(0) 



h = pC £i fturVei 

= pC e -VturC p Vf 



(10) 
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if the specific heat c p is assumed to be locally con- 
stant. The constant C £i is of order unity; it has been 
calibrated by numerical tests to provide the correct 
flame speed (see Niemeyer 1995 for details). Equa- 
tion ( p~0| ) shows that a temperature jump spreads over 
the distance A by turbulent diffusion during a char- 
acteristic time Tdiff w A 2 /j/ tU r ~ A/v(A). Since the 
burning timescale of a turbulent flame on the scale A 
that moves at a speed of u(A) ss v(A) is equally given 
by A/v(A) we see that the proposed flame model is 
self-consistent in the turbulent burning regime: the 
boundary of the high temperature region moves with 
the speed of a turbulent flame front. 

3. Simulations of central and off-center igni- 
tions 

The simulations involved an Eulerian PPM-based 
code to solve the two-dimensional hydrodynamical 
equations. Specifically, we employed the program 
PROMETHEUS ( [Fryxell et al. 1989| ). Its equation 
of state included contributions of ideal baryon and 
photon gases, and of relativistically degenerate elec- 
trons. The gravitational potential was computed in 
monopole approximation, while adding the resulting 
acceleration in a time symmetric fashion provided sec- 
ond order accuracy. 

In order to account for turbulent flame propaga- 
tion, the pro gram was modified with the S G-model 



described in Niemeyer fc Hillebrandt 1995. All cal- 



culations were based on a two-dimensional stationary 
grid with 256 x 64 zones in spherical (r, -d) coordi- 
nates. Assuming rotational and equatorial symmetry, 
the boundary conditions were chosen to be reflecting 
everywhere except at the outer radial edge, where out- 
flow was allowed. Owing to complications with unnat- 
urally high production of SG-turbulence at reflecting 
walls, both the SG-model and thermonuclear burning 
were inhibited in close vicinity of the grid boundaries. 
The zones were placed at equidistant radial and an- 
gular intervals with the exception of the 40 outermost 
ones whose radial separation was defined to increase 
at a rate of 10 % each. A Cray Y-MP computer at the 
Rechenzentrum Garching was used for all simulations. 

The initial model represented a M c h a n white dwarf 
with a central density p c = 2.8 x 10 9 g cm -3 and 
a central temperature T c = 7 x 10 8 K. In addition 
to initial temperature and density profiles a starting 
value for the turbulent kinetic SG-energy q had to be 
assigned. Earlier experiments showed that the calcu- 



lations are rather insensitive to the choice of go since 
the SG-model adjusts to the strength of local shear on 
very short timescales ( Niemeyer fc Hillebrandt 1995| ). 
A reasonable value is provided by the typical convec- 
tive velocities v conv « 5 x 10 5 cm/s on the grid scale 
A « 10 6 cm during the pre-ignition phase; this corre- 
sponds to q 10 11 erg/g. 

The flame front was initiated by instantly chang- 
ing the composition of a certain number of zones to 
nickel while increasing their temperature to several 
billion degrees, equivalent to raising their internal en- 
ergy by the binding energy e nuc m 7 x 10 17 erg/g. 
Immediately after the simulation has started pres- 
sure equalizes across the interface of burned and un- 
burned material by emitting weak compression and 
expansion waves in opposite directions. Meanwhile, 
turbulent transport increases the temperature in un- 
burned matter surrounding the burned regions and 
lowers their burning times. After a brief period of 
burning in the intermediate regime (no flamelets on 
small scales, transport dominated by turbulence on 
all scales) a conductive flame front forms after ap- 
proximately 10~ 3 s. 

3.1. Centrally ignited deflagrations 

In the central ignition (CI) model, a spherical re- 
gion around the star's center was ignited at t = 
s. Its surface was perturbed by a periodic pattern 
with a wavelength of 7r/4 and an amplitude of 25 km; 
this corresponded to a maximum radial distance of 
the initial flame front from the star's center of 150 
km. The development of the burning front in this 
and all other models is visualized by snapshots taken 
at 0.1, 0.8, and 1.4 s. In all figures, the velocity field is 
represented by randomly selected vector arrows at the 
respective points, superimposed with contour lines for 
the density and shaded areas to indicate the energy 
generation rate. 

Shortly after ignition, the flame propagates in a di- 
rection normal to its surface and forms a pa ttern that 
resem ble s the cellular structure discu ssed in Khokhlov 
1995| and [Niemeyer fc Wooslcy 1996] The cusps point- 



ing toward the burned matter are clearly visible in 
fig. (||). Large vortices have developed, causing the 
onset of a turbulent cascade that corresponds to the 
production of turbulent SG-energy in this model (note 
that SG-turbulence is only produced in non-burning 
zones). At t « 0.1 s, the turbulent regime dominates 
over laminar propagation almost everywhere within 
the front. 
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Fig. 1. — Scaling of conductive and turbulent length 
scales with the burning time. The shaded region 
marks the dominant burning length scale. 
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Fig. 2. — Core region in the Cl-model shortly after ig- 
nition. Density contours are separated by Sp = 5 • 10 7 
g cm -3 , the maximum energy generation rate (shaded 
regions) is S max = 10 19 erg g^ 1 s -1 . 



As the density encountered by the flame brush de- 
creases and effective gravity rises, buoyancy effects 
begin to govern the front evolution. The cellular 
structure is disrupted by the formation of hot bub- 
bles floating into the cold material. 

At t ss 0.8 s, displayed in fig. (|^), the RT- instability 
of the burning front has developed to its maximum ex- 
tent. The turbulent flame speed is now u(A) w2x 10 7 
cm/s at the points of maximum turbulent subgrid 
energy. Four major bubbles can be identified, sep- 
arated by thin streams of unburned C+O-material. 



As Khokhlov f994 points out, this particular pattern 
is an artifact of the two-dimensional simulation, as 
opposed to an increased fragmentation of the front 
in three-dimensional calculations. Since the mate- 
rial velocity on large scales is governed by buoyancy 
and the small scales behave as if fragmented in three- 
dimensional turbulence by virtue of the SG-model, 
the consequences of this effect may not be essential. 
However, as more efficient fragmentation also leads 
to an increased production of turbulence it is possible 
that it indirectly alters the strength of the explosion. 
It is still necessary to confirm this assertion by means 
of three-dimensional simulations. 

After the star has expanded and the density of the 
burning region has fallen to some 10 7 g cm~ 3 , the 
flame is extinguished at low densities. In our code, 
this occurs when /bum ~ A. At this stage, incom- 
plete nuclear burning presumably produces interme- 
diate mass elements that can be observed in SN la 
spectra. The simple network employed in this work is 
unable to predict the final composition of the burning 
products at densities below 10 7 g cm~ 3 . Therefore, 
the exact results of this stage of our simulations have 
to be interpreted with some care, since they can only 
qualitatively represent the physically realistic behav- 
ior of a dying flame. However, the energetical aspects 
of the explosion are widely unaffected by this short- 
coming and may therefore be taken at face value. It 
is also possible to give an estimate of the amount of 
intermediate mass elements produced in our models 
by keeping track of the amount of burned matter as 
a function of burning density (see section ^) . 

3.2. Off-center ignition at one point 

GW's arguments suggest that flame ignition may 
occur on the surface of buoyant bubbles that have al- 
ready risen to radial distances of approximately 200 
km . According to their estimates, the bubbles have 
diameters of some 10 km and velocities around 100 
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km/s at the time of flame formation. In order to 
study the effects of off-center ignition we carried out 
a number of simulations where the flame was ignited 
at only one point located at a radius of 200 km. As 
only one zone was initially burned in these one bubble 
(lB)-models, the radius of the ignited region roughly 
equaled 0.5A « 5 km. Some trial calculations in- 
volved giving the ignited zone an initial velocity of 
100 km/s, but after very short times the velocity field 
relaxed to the same solution as in the case of an ini- 
tially resting ignition zone. It should be kept in mind 
in the ensuing discussion that the assumption of rota- 
tional and equatorial symmetry enters strongly in this 
simulation. Despite the fact that the actual structure 
of the burning front represents two toroidal buoyant 
rings we shall keep the term "bubble" in the following 
description of the front evolution. 

The bubble quickly attains a drop-like shape and 
instantly begins to float with some I0 7 cm/s (fig. |). 
The energy generation rate is weak compared with 
the Cl-model at the same time, resulting from the 
lack of strong shear flows. This can be explained by 
the larger amount of potential energy of hot material 
in the Cl-model at the time of ignition. It is quickly 
transformed into kinetic energy as the RT-structures 
begin to float, giving rise to the production of turbu- 
lent SG-energy. The floating volume in the IB-model, 
on the other hand, is smaller in the beginning and 
therefore unable to provide a large amount of turbu- 
lent kinetic energy to accelerate the burning front. 

Having reached a maximum radial extent of ap- 
proximately 450 km and a rising velocity of more than 
10 s cm/s, the bubble begins to spread out in the an- 
gular dimension, fig. (||) shows the burning bubble 
in the fully developed "mushroom stage" : a nearly 
spherical head with a diameter of roughly 500 km is 
attached to a thin pole that connects it to the center 
of the expanding star. The most prominent feature 
of this configuration is the large volume of burning 
material at low densities compared to the part of the 
front that burns at higher density. It may help to 
account for the production of intermediate mass ele- 
ments necessary to reproduce SN la spectra. A more 
detailed discussion of this aspect will follow in section 




3.3. Ignition at five points 

As only little is known about the circumstances of 
flame ignition it is also possible that ignition occurs 
at many points scattered throughout the white dwarf 
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Fig. 3.- Sp = 3 • I0 7 g cm" 3 , S max = 1.3 • I0 19 erg 
g" 1 s- 1 . 
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Fig. 4. — IB-model with dp = 6 ■ I0 7 g cm 3 and 
Smax^-IO^ergg- 1 s~\ 
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core region almost simultaneously. The multiple igni- 
tion model (specifically, ignition at five points (5B)) 
can be viewed as a compromise between the instan- 
taneous incineration of the core region at t = and 
ignition at only one point. This way, some desirable 
properties of both the Cl-model (strong velocity gra- 
dients triggering a powerful turbulent cascade) and 
the IB-model (burning at lower density, increasing 
the synthesis of intermediate mass elements) may be 
combined. 

The model presented here differs from the standard 
models described for the CI and IB-cases with respect 
to the turbulent diffusion coefficient C £i (section |^) . 
While C Ci — 2 was used in all other simulations of tur- 
bulence dominated flames this calculation was carried 
out with C Ci = 1. It therefore models a sligthly more 
slowly propagating flame front, while the burning rate 
remained unchanged. In order to assess the influence 
of the propagation and burning rates on the explo- 
sion a second simulation with identical initial condi- 
tions but different calibration parameters (C ei = 2 
and B = 0.511) was carried out; in section (||), these 
models are referred to as the "slow" and "fast" SB- 
models, respectively. 

The initial ignition pattern chosen for the SB- 
model can be seen in fig. (||) . Three zones were incin- 
erated at r w 100 km while two more bubbles were 
triggered at r ~ 200 km. This choice follows from 
the intuitive notion that ignition is more likely at 
higher temperatures closer to the center. Initially, 
all five bubbles evolve independently, floating in the 
same way as described for the IB-model. After a short 
time, the bubbles begin to interact. As a consequence 
of their lower buoyant speed the inner three bubbles 
merge first, later joined by by the upper two. 

It is interesting to compare fig. ([?]) with figs, 
and (^). While the maximum material velocity 
(shown at the bottom of the plots) of the 5B-modcl 
equals that of the central ignition calculation, its cen- 
tral density is considerably lower. Furthermore, most 
of the burning material is further off center and there- 
fore burns at a lower average density, similar to the 
IB-model. However, the energy generation rate in 
the 5B-model is higher than in the latter, indicating 
stronger turbulent motions. The apparent similar- 
ity of front evolution in both off-center models and 
their difference from the Cl-simulation gives evidence 
that under most natural initial conditions, the bulk 
volume of burning material can be expected to float 
upward. The fact that this is not observed in central 
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Fig. 5.— dp = 10 7 g cm- 3 , S max = 10 19 erg g" 1 s" 1 . 
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Fig. 6. — "Slow" 5B-modcl; dp = 6 • 10 7 g cm 3 and 
£ ma x = 6-10 18 ergg- 1 s~\ 
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ignition simulations is probably related to the sym- 
metry conditions that are usually imposed at the grid 
boundaries. 

At the time of flame extinction, the high extent 
of fragmentation compared with the other models be- 
comes visible. The flame covers a substantial part 
of the picture and still burns strongly in regions of 
higher density. This fact explains the higher release 
of energy in the 5B-explosion compared with the other 
models (section^). 

4. Global energetic parameters 

The most important parameter to determine the fi- 
nal fate of the white dwarf is its total energy. So far, 
all multidimensional simulations of the problem have 
failed to release the required ;> 0.5 M Q in the first 
deflagration phase needed to account for SN la ob- 
servations, therefore implying pulsational detonations 
if Chandrasekhar mass models were to explain these 



.Time 



events (Khokhlov 1995, Arnett & Livne 1994a, Arnett 
Livne 1994b). The simulations presented in this 
work provide a possible explosion mechanism with- 
out involving detonations, but the explosions are still 
too weak to explain most Type la supernovae. This 
can be seen in fig. (|8|), showing the total energy as a 
function of time. All models start at a value of ap- 
proximately —5 x 10 50 erg, consisting mostly of neg- 
ative gravitational energy and positive Fermi energy 
of the degenerate electron gas. Note that the CI- 
model (solid line) has a slightly higher total energy 
at t = which is related to the larger initially incin- 
erated mass. 

The only model that remains bound after burning 
has ceased is the IB-simulation. It burns less than 
0.3 M in the first deflagration phase and is thus the 
only candidate for a pulsational detonation. Accord- 
ing to the description in section (B.2), the front fails to 
ignite a sufficiently large volume because of its lim- 
ited volume and the lack of fragmentation. While 
the latter aspect may change in a three dimensional 
treatment, it must be stressed that the simulation al- 
ready assumes rotational and equatorial symmetry. 
Three-dimensional simulations of the full star would 
therefore show a much smaller volume of the burning 
bubble and thus predict the generation of consider- 
ably less energy. Hence, it is unlikely that one-point 
ignitions yield enough energy for a healthy explosion 
on the first try unless the burning front spreads sig- 
nificantly faster than modeled. 



e+01 !,p,= 4,93e+0§F a/cm" 
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Fig. 7.— Sp = 3 • 10 7 g cm- 3 , S max = 10 19 erg g" 1 
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Fig. 8. — Total energy of the different models as a 
function of time ( — CI, - - IB, - • - slow 5B, - • • • fast 
5B). 
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One may, of course, speculate about the subse- 
quent evolution of the IB-model. Some burning mate- 
rial, surrounded by large amounts of unburned carbon 
and oxygen, probably remains in the flamelet regime 
until the star recontracts. Re-ignition could then oc- 
cur on a highly convoluted and aspherical surface 
that results from folding and compression of the thin 
stream of burning matter reaching toward the center. 
Whether this process is able to form a detonation, or 
else produces more buoyant bubbles with turbulent 
burning surfaces cannot be answered at this point. 
In order to obtain information about the re-ignition 
phase the multidimensional evolution has to be fol- 
lowed for several seconds after the initial deflagration, 
which is computationally unfeasible at present. 

The Cl-explosion terminates burning after produc- 
ing 0.59 M Q of nuclear fuel, which represents the en- 
ergetical threshold of an observable, but weak, SN la 
event. Despite its higher mass of initially burned ma- 
terial and the relatively large surface of the burning 
front it releases less energy than the 5B-models. 

Finally, both fast and slow 5B-simulations appear 
to represent possible candidates for some Type la su- 
pcrnovac. Independent of the variations in burning 
and propagation speed mentioned in section (3.3), be- 
tween 0.65 M and 0.67 M Q matter are burned, lib- 
erating approximately 9 x 10 50 erg of nuclear energy 
(this number is sligthly overestimated due to our as- 
sumption of NSE even at low burning densities, see 
the next paragraph). 

To conclude the discussion of the performed nu- 
merical experiments some implications for nucleosyn- 
thesis in SN la events shall be pointed out. On ac- 
count of the simplified nuclear network that implies 
NSE for the burning products, only indirect state- 
ments about the expected final composition are pos- 
sible. For this reason, the amount of burned carbon 
and oxygen as a function of burning density is plot- 
ted in fig. for all models. Of particular interest is 
the region below 10 7 g cm -3 where intermediate mass 
elements begin to be synthesized. The solid curve, 
representing the Cl-model, shows a sharp peak at ap- 
proximately 4 x 10 8 g cm' 3 followed by a rapid decline 
toward lower densities. Very little burning occurs in 
the intermediate mass region. On the other hand, all 
three off-center simulations burn significantly larger 
amounts of matter below 10 7 g cm -3 . Most surprising 
is the ratio of this value to the total amount of burned 
matter for the IB-model: it reaches 20 percent, com- 
pared to only four percent for the Cl-model. Addi- 



tional information about the expected nucleosynthe- 
sis can be gained from the region above w 1.5 x 10 9 
g cm~ 3 . There, the high electron density leads to 
excess neutronization compared to solar composition 
(Woosley 1990). As expected, the IB-model burns 



the smallest amount of C+O-matter at high densi- 
ties. However, both 5B-models burn even more than 
the Cl-simulation which might be explained by the 
rapid initial expansion of the latter during the first 
10 -2 s because of its higher mass of initially ignited 
material. 

The total mass of burned material near the end of 
the simulations is summarized in table (^), together 
with the parameters of equations (^) and ( |l0| ) used 
for the different models. 

5. Conclusions 

It can reasonably be assumed that flame ignition 
in M c han models for white dwarf explosions occurs 
under convective conditions that make ignition at the 
precise center unlikely. Based upon a simple model for 
the ignition process (GW) , we expect that in general, 
ignition will occur at one or many points within the 
core. Three models were studied in two-dimensional 
simulations: ignition at one (IB) and five (5B) points 
slightly off-center, in addition to the standard central 
ignition (CI) scenario. 

The overall appearance of both off-center models 
is strikingly different from the Cl-case. As the bub- 
ble (or its merging components) rises, it leaves behind 
a trail of burning material while the uppermost part 
grows rapidly. At the maximum of global energy gen- 
eration the shape of the burning front resembles a 
mushroom whose wide head is supported by a thin 
pole that reaches from the center to about a third of 
the star's radius. Various processes are responsible 
for this structure: first, the turbulent flame speed is 
much less at the bottom because of the smaller ve- 
locity gradients and shear stress. Therefore, angular 
spreading at the pole happens at a smaller rate than 
in the stirred upper region. Equally important is the 
flow pattern around and within the bubble: displace- 
ment of cold material by the rising volume creates 
large vortices that diverge at the top and converge at 
the bottom of the bubble, thereby focusing the pole 
and expanding the top. They also cause the high 
burning velocities in the upper part by creating high 
velocity gradients in the mixing region. Formation 
of the macroscopic vortices therefore also determines 
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the radius (and, consequently, density) at which the 
bubble begins to spread out. In the simulations, this 
happens at approximately one third to one half of 
the star's radius. It must be kept in mind that the 
details of this development, in particular the forma- 
tion of large scale vortices, are probably different in 
three-dimensional reality. 

The RT-fingers of the Cl-model display similar 
characteristics on smaller scales, but never develop 
the same elongated mushroom shape. This is at least 
partly due to the fact that initial perturbations first 
have to overcome the linear RT growth stage before 
becoming fully buoyant. By the time the nonlinear 
evolution sets in, the dynamics of the largest scales 
is already strongly influenced by the expansion of 
the star. From this point of view, off-center igni- 
tion models simply skip the linear instability regime 
and consequently gain important time to rise and be- 
come turbulent. The most important consequence of 
this behavior is a higher total burning rate at early 
times than predicted by standard deflagration scenar- 
ios. We can expect this effect to be even more pro- 
nounced in three dimensions and with a large number 
of initial ignition points. If ignition occurcd almost 
simultaneously throughout the central region, but lo- 
cated on bubble surfaces with small radii compared 
to that of the ignition region, the bubbles would reach 
their maximum burning rate at roughly the same 
time. The amount of matter that has been burned 
by then, which in turn determines the degree of ex- 
pansion of the star at the time rapid burning sets in, 
will have to be calculated by three-dimensional simu- 
lations with high resolution. 

A rough estimate of the average bubble velocity in 
the limit of many interacting bubbles in three dimen- 



sions can be obtained from experiments (Read 1984 



snider & Andrews 1994) and simulations (Youngs 
198;4|, |Li 1996| ) of RT mixing fronts in the fully non- 
linear regime. It was shown that the front boundary 
propagates into the unmixed fluid with the velocity 



Wfront = 2 ag e s t , 



(11) 



independent of the initial conditions. Here, t is the 
time since the beginning of the nonlinear RT stage 
and a « 0.07 is a constant. This behavior was 
successfully explained by statistical models pioneered 



by Sharp and Wheeler (Sharp 1984) and refined by 
Slimm fc Li 1988| , |Glimm fe Sharp 199C| and [Zhang 
19SjO|. The Sharp- Wheeler model assumes that inter- 
acting bubbles merge to form bigger ones while small 



bubbles are being washed downstream, thus increas- 
ing the average bubble radius (r) at the front bound- 
ary. Together with the rise velocity of single bubbles 
(0), this process can account for a constant increase 
of Uf r0 nt- It is easily seen that (^) and ( |TT] ) agree if 
(r) grows linearly with the distance of the bubbles 
from the inner front boundary, which is in our case 
equal to the size of the burned out region (this is 
also clear from dimensional analysis). Furthermore, 
our bubbles increase their volume by rapid burning 
and by gas expansion, providing additional growth as 
they rise outward. Hence, equation ( |ll|) is a reason- 
able approximation even if bubble merging is not the 
dominant dynamical effect. Using typical numbers 
for white dwarf explosions, <? c ff ~ 10 9 cm s~ 2 and 
i»1...2s,we obtain front velocities of the order of 
2 x 10 8 cm s" 1 . According to one dimensional stud- 



ies ( Woosley 1990 ) with similar velocities, this front 
velocity is sufficient to produce a healthy explosion. 

As a second result of off-center ignition, consider- 
ably more matter is burned in a lower density envi- 
ronment than in the central ignition simulation (ta- 
ble [l]). This helps to satisfy two major constraints 
on potential SN la explosion mechanisms: only lit- 
tle nuclear burning at high density is permitted to 
limit the overproduction of neutron-rich isotops, and 
a large amount of intermediate mass nuclei produced 
at p < 10 7 g cm" 3 is required to explain the ob- 



served spectral lines (Woosley 1990). Comparing the 
respective ratios of both masses and the total burned 
matter, the IB-model is clearly the best candidate in 
this particular respect. While the 5B-models produce 
large amounts of intermediate mass isotops they also 
burn more high density material than the Cl-model. 
In the overall comparison of the expected nucleosyn- 
thesis, however, all off-center models are closer to the 
expectations fr om SN la obse rvations and solar abun- 
dance studies ( Woosley 1990 ) than the Cl-model. 

Looking at the final values for the total energy the 
IB-model fails to be a Type la supernova (neglecting 
the possibility of a pulsational detonation). Despite 
the symmetry assumptions that imply identical en- 
ergy generation in both hemispheres, too little nuclear 
energy is released to unbind the star. We can there- 
fore conclude that in order to produce a successful 
explosion in the first deflagration phase, flame igni- 
tion must occur at a number of points within a short 
period of time (less than w 10 _1 s when the star be- 
gins to expand) , where the term "points" refers to re- 
gions with a diameter equal to or less than 10 km (the 
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size of our grid zones). Both 5B-simulations and the 
CI- model become unbound during their explosions. 
However, the Cl-case would only represent a weak SN 
la since it barely produces the minimum value of ap- 
proximately 0.6 M Q (table |l|) of nuclear burning prod- 
ucts necessary to account for light curves and kinetic 
energy of the ejecta (e.g., |Woosley 1990| ). The SB- 
calculations, on the other hand, are in most respects 
possible candidates for some SN la events, although 
the final amount of produced 56 Ni is still insufficient 
to account for the light curves of the full class of ex- 
plosions. We can thus summarize that deflagration 
models of M c h an white dwarfs with a small number 
of burning bubbles are unlikely to produce powerful 
SN la's unless our simulations fail to cover the rele- 
vant physics. In addition to the expected changes in 
three dimensions described above, examples for po- 
tentially important, but neglected effects include the 
additional production of turbulence by burning itself 
(Nicmeyer & Woosley 1996) or from differential rota- 
tion of the star. 

Of course, all results of supernova simulations have 
to be placed in the context of their numerical tech- 
nique. In our case, two numerical tools, the SG 
model for unresolved turbulence and the flame prop- 
agation algorithm, deserve closer inspection as they 
are both new and crucial ingredients in our work. 
The SG model and some of its numerical properties 
have been discussed in (Nicmeyer & Hillebrandt 1995 
and references therein). Within its limits (no pro- 
duction of turbulence by thermal expansion within 
the flame brush, artificial closure assumptions) it is 
capable of representing the unresolved turbulence in 
a way that we consider realistic. Dramatic changes 
in the outcome of the simulations by changing the 
SG para meters within reasonable limits cannot be ex- 
pected ( Niemeyer fc Hillebrandt 1995 ). The same 
is true for the fla me algorithm, whose performance 
has been tested in Niemeyer 1995 . However, more ef- 
fort is necessary to study the problem of numerically 
representing turbulent flame fronts in finite difference 
schemes like PPM. So far, the turbulent flame brush 
has either been treated like a thin, laminar flame sur- 
face that moves with the turbulent burning speed 
dKhokhlov 1995|, [Niemeyer fc Hillebrandt 1995[), or 



although it is unclear in which direction. 

The possible transition to a delayed or, in case of 
the IB model, pulsational detonation has not been 
considered in this work but remains, in principle, a 
viable option for SN la models. However, all of our 
results show subsonic effective burning velocities at all 
times during the turbulent deflagration phase, making 
the delayed formation of a detonation very unlikely. 
Again, the situation may turn out differently if one 
consideres an ensemble of many bubbles that start 
interacting almost simultaneously. 

The sensitivity of our results to the initial condi- 
tions opens a new dimension of parameter space that 
is worth exploring. One can speculate that some of 
the variations noticed in SN la observations can be 
accounted for by this effect. However, if ignition oc- 
curs at sufficiently many points throughout the core 
region, the outcome is again expected to be statisti- 
cally homogeneous. Detailed studies of the ignition 
process will be necessary to answer this question. 
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has been smeared out over a large region by turbu- 
lent diffusion (this work) . Both methods are probably 
wrong, at least to some extent. It cannot be excluded 
that increasing the numerical resolution and our abil- 
ity to model turbulent flames will change our results, 
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Table 1 

subgrid parameters and burned masses after t « 1.4 s for all models. 



Model 


B 




M b [M ] 


M b (p ;> 1.5 • 10 9 ) [M ] 


M b ( P ^10 7 ) [Mq] 


CI 


0.25 


2 


0.59 


3.5 ■ 10~ 2 


2.5 • 10~ 2 


IB 


0.25 


2 


0.27 


1.1 ■ 10~ 2 


5.3 • 10~ 2 


5B, slow 


0.25 


1 


0.65 


4.5 • 10~ 2 


6.9 • 10~ 2 


5B, fast 


0.511 


2 


0.67 


4.8- 10~ 2 


8.3 • 10~ 2 




Fig. 9. — Mass of burned matter as a function of the 
burning density at t w 1.4 s ( — CI, - - IB, - • - slow 
5B, fast 5B). 



13 



REFERENCES 

Arnctt, W. D., & Livnc, E. 1994a, ApJ, 427, 315 

Arnett, W. D., & Livne, E. 1994a, ApJ, 427, 330 

Bcnz, W. 1996, in Thermonuclear Supernovae, eds. 
R. Canal & P. Ruiz LaPuentc, in press. 

Fryxell, B.A., Miiller, E., & Arnctt, W.D. 1989, MPA 
Preprint Scries, 449 

Garcia-Scnz, D., Bravo, E., & Woosley, S. E. 1996, in 
Thermonuclear Supernovae, eds. R. Canal and P. 
Ruiz-LaPuente, in press 

Garcia-Scnz, D., & Woosley, S. E. 1995 (GW), ApJ, 
454, 895 

Glimm, J., & Li, X.L. 1988, Phys. Fluids, 31, 2077 

Glimm, J., & Sharp, D.H. 1990, Phys. Rev. Lett., 64, 
2137 

Hoflich, P., & Khokhlov, A. 1996, ApJ, in press. 

Iben, I., Jr. 1982, ApJ, 253, 248 

Iben, I., Jr., & Tutukov, A. V. 1991, ApJ, 370, 615. 

Kcnyon, S. J., Livio, M., Mikolajewska, J., & Tout, 
C. A. 1993, ApJ, 407, L81 

Khokhlov, A.M. 1991a, A&A, 245, 114 

Khokhlov, A.M. 1991b, A&A, 245, L25 

Khokhlov, A.M. 1991c, A&A, 246, 383 

Khokhlov, A.M. 1993, ApJ, 419, L77 

Khokhlov, A.M. 1994, ApJ, 424, L115 

Khokhlov, A.M. 1995, ApJ, 449, 695 

Landau, L.D., & Lifshitz, E.M. 1991, Lchrbuch der 
Thcoretischen Physik VI: Hydrodynamik (Berlin: 
Akademie-Verlag) 

Layzer, D. 1955, ApJ, 122, 1 

Li, X.L. 1996, Phys. Fluids, 8, 336 

Limongi, M., & Tornambc, A. 1991, ApJ, 371, 317. 

Livne, E. 1993, ApJ, 406, L17 

Livne, E., & Arnett, W.D. 1996, in Thermonuclear 
Supernovae, eds. R. Canal & P. Ruiz LaPuente, in 
press 



Miiller, E. 1994, Fundamentals of Gasdynamical Sim- 
ulations, MPA Preprint Series, 780 

Miiller, E., & Arnett, W.D. 1986, ApJ, 307, 619 

Munari, IT., & Renzini, A. 1992, ApJ, 397, L87 

Niemeyer, J.C. 1995, Ph.D. Thesis, MPA Preprint Se- 
ries, 911 

Niemeyer, J.C, & Hillebrandt, W. 1995, ApJ, 452, 
769 

Niemeyer, J.C, & Woosley, S.E. 1996, subm. to ApJ 

Nomoto, K., Thielemann, F.-K., & Yokoi, K. 1984, 
ApJ, 286, 644 

Rappaport, S., Di Stefano, R., & Smith, J. D. 1994, 
ApJ, 426, 692 

Read, K.I. 1984, Physica D, 12, 45 

Sharp, D.H. 1984, Physica D, 12, 3 

Snider, D.M., & Andrews, M.J. 1994, Phys. Fluids, 
6, 3324 

Timmcs, F.X., & Woosley, S.E. 1992, ApJ, 396, 649 

Woosley, S.E. 1990, in Supernovae, ed. A.G. Petschck 
(Berlin: Springer- Verlag) , p. 182 

Woosley, S. E., and Eastman, R. G. 1995, Proceedings 
of Menorca School of Astrophysics, eds. E. Bravo, 
R. Canal, J. Ibanez, and J. Isern, Societat Catalana 
de Fisica, p. 105 

Woosley, S. E., & Weaver, T. A. 1986, in Radiation 
Hydrodynamics in Stars and Compact Objects, 
eds. D. Mihalas & K.-H. A. Winkler, (Springer Ver- 
lag:Berlin), p. 91 

Woosley, S. E., & Weaver, T. A. 1994a, ApJ, 423, 
371. 

Youngs, D.L. 1984, Physica D, 12, 32 
Zhang, Q. 1990, Phys. Lett. A, 151, 18 



This 2-column preprint was prepared with the AAS IATj^X 
macros v4.0. 



14 



